{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Integration and Differentiation\n",
    "\n",
    "This notebook illustrates how to perform numerical integration and differentiation. \n",
    "\n",
    "There are several packages for doing this. Here, the focus is on [QuadGK](https://github.com/JuliaMath/QuadGK.jl) and [ForwardDiff](https://github.com/JuliaDiff/ForwardDiff.jl)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "printyellow (generic function with 1 method)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Dates, QuadGK, ForwardDiff\n",
    "\n",
    "include(\"printmat.jl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "\n",
    "#pyplot(size=(600,400))\n",
    "gr(size=(480,320))\n",
    "default(fmt = :svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Integration\n",
    "\n",
    "As a simple illustration, the next cells plot and integrate the $N(0,\\sigma)$ pdf."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Pdf of the N() Distribution \n",
    "\n",
    "Notice: the function takes $\\mu$ and $\\sigma$ (not $\\sigma^2$) as inputs, similar to the Distributions package. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ϕNS (generic function with 3 methods)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function ϕNS(x,μ=0,σ=1)       #pdf of N(μ,σ), defaults to N(0,1)\n",
    "    z   = (x - μ)/σ\n",
    "    pdf = exp(-0.5*z^2)/(sqrt(2*pi)*σ)\n",
    "    return pdf\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"480\" height=\"320\" viewBox=\"0 0 1920 1280\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip3300\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"1920\" height=\"1280\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip3300)\" d=\"\n",
       "M0 1280 L1920 1280 L1920 0 L0 0  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip3301\">\n",
       "    <rect x=\"384\" y=\"0\" width=\"1345\" height=\"1280\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip3300)\" d=\"\n",
       "M153.898 1105.62 L1872.76 1105.62 L1872.76 121.675 L153.898 121.675  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip3302\">\n",
       "    <rect x=\"153\" y=\"121\" width=\"1720\" height=\"985\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip3302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  202.545,1105.62 202.545,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  472.806,1105.62 472.806,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  743.066,1105.62 743.066,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1013.33,1105.62 1013.33,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1283.59,1105.62 1283.59,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1553.85,1105.62 1553.85,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1824.11,1105.62 1824.11,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  153.898,1077.77 1872.76,1077.77 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  153.898,845.094 1872.76,845.094 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  153.898,612.416 1872.76,612.416 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  153.898,379.739 1872.76,379.739 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  153.898,147.062 1872.76,147.062 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.898,1105.62 1872.76,1105.62 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.898,1105.62 153.898,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  202.545,1105.62 202.545,1090.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  472.806,1105.62 472.806,1090.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  743.066,1105.62 743.066,1090.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1013.33,1105.62 1013.33,1090.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1283.59,1105.62 1283.59,1090.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1553.85,1105.62 1553.85,1090.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1824.11,1105.62 1824.11,1090.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.898,1077.77 179.681,1077.77 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.898,845.094 179.681,845.094 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.898,612.416 179.681,612.416 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.898,379.739 179.681,379.739 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.898,147.062 179.681,147.062 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 202.545, 1157.22)\" x=\"202.545\" y=\"1157.22\">-3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 472.806, 1157.22)\" x=\"472.806\" y=\"1157.22\">-2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 743.066, 1157.22)\" x=\"743.066\" y=\"1157.22\">-1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1013.33, 1157.22)\" x=\"1013.33\" y=\"1157.22\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1283.59, 1157.22)\" x=\"1283.59\" y=\"1157.22\">1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1553.85, 1157.22)\" x=\"1553.85\" y=\"1157.22\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1824.11, 1157.22)\" x=\"1824.11\" y=\"1157.22\">3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 134.698, 1095.27)\" x=\"134.698\" y=\"1095.27\">0.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 134.698, 862.594)\" x=\"134.698\" y=\"862.594\">0.1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 134.698, 629.916)\" x=\"134.698\" y=\"629.916\">0.2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 134.698, 397.239)\" x=\"134.698\" y=\"397.239\">0.3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 134.698, 164.562)\" x=\"134.698\" y=\"164.562\">0.4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:84px; text-anchor:middle;\" transform=\"rotate(0, 1013.33, 73.2)\" x=\"1013.33\" y=\"73.2\">pdf of N(0,1)</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1013.33, 1234.68)\" x=\"1013.33\" y=\"1234.68\">x</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip3302)\" style=\"stroke:#ff0000; stroke-width:8; stroke-opacity:1; fill:none\" points=\"\n",
       "  202.545,1067.46 229.571,1063.92 256.597,1059.35 283.623,1053.52 310.649,1046.17 337.675,1036.99 364.701,1025.66 391.727,1011.86 418.753,995.23 445.78,975.431 \n",
       "  472.806,952.146 499.832,925.098 526.858,894.072 553.884,858.94 580.91,819.683 607.936,776.413 634.962,729.389 661.988,679.035 689.014,625.944 716.04,570.878 \n",
       "  743.066,514.76 770.092,458.651 797.118,403.724 824.144,351.227 851.171,302.433 878.197,258.595 905.223,220.89 932.249,190.368 959.275,167.903 986.301,154.152 \n",
       "  1013.33,149.523 1040.35,154.152 1067.38,167.903 1094.41,190.368 1121.43,220.89 1148.46,258.595 1175.48,302.433 1202.51,351.227 1229.54,403.724 1256.56,458.651 \n",
       "  1283.59,514.76 1310.61,570.878 1337.64,625.944 1364.67,679.035 1391.69,729.389 1418.72,776.413 1445.74,819.683 1472.77,858.94 1499.8,894.072 1526.82,925.098 \n",
       "  1553.85,952.146 1580.87,975.431 1607.9,995.23 1634.93,1011.86 1661.95,1025.66 1688.98,1036.99 1716,1046.17 1743.03,1053.52 1770.06,1059.35 1797.08,1063.92 \n",
       "  1824.11,1067.46 \n",
       "  \"/>\n",
       "<path clip-path=\"url(#clip3302)\" d=\"\n",
       "M202.545 1067.46 L229.571 1063.92 L256.597 1059.35 L283.623 1053.52 L310.649 1046.17 L337.675 1036.99 L364.701 1025.66 L391.727 1011.86 L418.753 995.23 L445.78 975.431 \n",
       "  L472.806 952.146 L499.832 925.098 L526.858 894.072 L553.884 858.94 L580.91 819.683 L607.936 776.413 L634.962 729.389 L661.988 679.035 L689.014 625.944 L716.04 570.878 \n",
       "  L743.066 514.76 L770.092 458.651 L797.118 403.724 L824.144 351.227 L851.171 302.433 L878.197 258.595 L905.223 220.89 L932.249 190.368 L959.275 167.903 L986.301 154.152 \n",
       "  L1013.33 149.523 L1040.35 154.152 L1067.38 167.903 L1094.41 190.368 L1121.43 220.89 L1148.46 258.595 L1175.48 302.433 L1202.51 351.227 L1229.54 403.724 L1256.56 458.651 \n",
       "  L1283.59 514.76 L1310.61 570.878 L1337.64 625.944 L1364.67 679.035 L1391.69 729.389 L1418.72 776.413 L1445.74 819.683 L1445.74 1077.77 L1418.72 1077.77 L1391.69 1077.77 \n",
       "  L1364.67 1077.77 L1337.64 1077.77 L1310.61 1077.77 L1283.59 1077.77 L1256.56 1077.77 L1229.54 1077.77 L1202.51 1077.77 L1175.48 1077.77 L1148.46 1077.77 L1121.43 1077.77 \n",
       "  L1094.41 1077.77 L1067.38 1077.77 L1040.35 1077.77 L1013.33 1077.77 L986.301 1077.77 L959.275 1077.77 L932.249 1077.77 L905.223 1077.77 L878.197 1077.77 L851.171 1077.77 \n",
       "  L824.144 1077.77 L797.118 1077.77 L770.092 1077.77 L743.066 1077.77 L716.04 1077.77 L689.014 1077.77 L661.988 1077.77 L634.962 1077.77 L607.936 1077.77 L580.91 1077.77 \n",
       "  L553.884 1077.77 L526.858 1077.77 L499.832 1077.77 L472.806 1077.77 L445.78 1077.77 L418.753 1077.77 L391.727 1077.77 L364.701 1077.77 L337.675 1077.77 L310.649 1077.77 \n",
       "  L283.623 1077.77 L256.597 1077.77 L229.571 1077.77 L202.545 1077.77  Z\n",
       "  \" fill=\"#ff0000\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip3302)\" style=\"stroke:#e26f46; stroke-width:8; stroke-opacity:1; fill:none\" points=\"\n",
       "  202.545,1067.46 229.571,1063.92 256.597,1059.35 283.623,1053.52 310.649,1046.17 337.675,1036.99 364.701,1025.66 391.727,1011.86 418.753,995.23 445.78,975.431 \n",
       "  472.806,952.146 499.832,925.098 526.858,894.072 553.884,858.94 580.91,819.683 607.936,776.413 634.962,729.389 661.988,679.035 689.014,625.944 716.04,570.878 \n",
       "  743.066,514.76 770.092,458.651 797.118,403.724 824.144,351.227 851.171,302.433 878.197,258.595 905.223,220.89 932.249,190.368 959.275,167.903 986.301,154.152 \n",
       "  1013.33,149.523 1040.35,154.152 1067.38,167.903 1094.41,190.368 1121.43,220.89 1148.46,258.595 1175.48,302.433 1202.51,351.227 1229.54,403.724 1256.56,458.651 \n",
       "  1283.59,514.76 1310.61,570.878 1337.64,625.944 1364.67,679.035 1391.69,729.389 1418.72,776.413 1445.74,819.683 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1486.28, 487.658)\" x=\"1486.28\" y=\"487.658\">the area covers</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1486.28, 539.498)\" x=\"1486.28\" y=\"539.498\">up to x=1.645</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x  = -3:0.1:3\n",
    "xb = x[x.<=1.645]\n",
    "\n",
    "plot( x,ϕNS.(x),\n",
    "      linecolor = :red,\n",
    "      linewidth = 2,\n",
    "      legend = nothing,\n",
    "      title = \"pdf of N(0,1)\",\n",
    "      xlabel = \"x\",\n",
    "      ylabel = \"\",\n",
    "      annotation = (1.75,0.25,text(\"the area covers\\nup to x=1.645\",8)) )\n",
    "\n",
    "plot!(xb,ϕNS.(xb),fillcolor=:red,linewidth=2,legend=nothing,fill=(0,:red))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating Prob(x<=1.645)\n",
    "\n",
    "The next cell calculates (by numeric integration)\n",
    "\n",
    "$\\int_{-\\infty}^{1.645}\\phi(x,0,\\sigma)dx$,\n",
    "\n",
    "where $\\phi(x,\\mu,\\sigma)$ is the pdf of an $N(\\mu,\\sigma)$ variable.\n",
    "\n",
    "\n",
    "The input to `quadgk` should be a function with only one argument. We do that by creating an anonymous function\n",
    "```\n",
    "x->fn(x,a),\n",
    "```\n",
    "assuming that `a` has a value already."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Pr(x<=1.64) according to N(0,1):     0.950\n",
      "\n",
      "Pr(x<=1.64) according to N(0,2):     0.795\n"
     ]
    }
   ],
   "source": [
    "cdf1, = quadgk(x->ϕNS(x),-Inf,1.645)       #N(0,1)\n",
    "printlnPs(\"\\nPr(x<=1.64) according to N(0,1):\", cdf1)\n",
    "\n",
    "cdf2, = quadgk(x->ϕNS(x,0,2),-Inf,1.645)   #N(0,σ=2)\n",
    "printlnPs(\"\\nPr(x<=1.64) according to N(0,2):\", cdf2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Numerical Derivatives\n",
    "\n",
    "Numerical derivatives can be calculated by a crude finite difference (see `NumDer()` below) or by the much more sophisticated routines in the `ForwardDiff` package.\n",
    "\n",
    "To do this calculation for many `x` values, we loop (either a `for` loop or list comprehension). Dot syntax, `ForwardDiff.derivative.()`, seems to work, but is not documented."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "NumDer (generic function with 1 method)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function NumDer(fun,b0,h)     #crude function for a centered numerical derivative\n",
    "    bminus = b0 .- h\n",
    "    bplus  = b0 .+ h\n",
    "    hh     = bplus - bminus\n",
    "    fplus  = fun(bplus)\n",
    "    fminus = fun(bminus)\n",
    "    D      = (fplus-fminus)/hh\n",
    "    return D\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "fn1 (generic function with 1 method)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function fn1(x,a)                 #a simple function, to be differentiated\n",
    "   return (x - 1.1)^2 - a\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The derivative at x=2 is (from two different methods): \n",
      "     1.800     1.800\n",
      "\n"
     ]
    }
   ],
   "source": [
    "x0 = 2\n",
    "\n",
    "dydx_A = NumDer(x->fn1(x,0.5),x0,0.01)           #differentiate fn1(x,0.5) at x = x0\n",
    "dydx_B = ForwardDiff.derivative(x->fn1(x,0.5),x0)\n",
    "\n",
    "println(\"The derivative at x=$x0 is (from two different methods): \")\n",
    "printmat([dydx_A dydx_B])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "now lets plot this\n"
     ]
    }
   ],
   "source": [
    "x = -3:6/99:6          #calculate the derivative at many points\n",
    "\n",
    "dydx_A = [NumDer(x->fn1(x,0.5),x[i],0.01) for i=1:length(x)]  #list comprehension as a quick loop\n",
    "\n",
    "dydx_B = [ForwardDiff.derivative(x->fn1(x,0.5),x[i]) for i=1:length(x)]\n",
    "\n",
    "println(\"now lets plot this\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"480\" height=\"320\" viewBox=\"0 0 1920 1280\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip3700\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"1920\" height=\"1280\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip3700)\" d=\"\n",
       "M0 1280 L1920 1280 L1920 0 L0 0  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip3701\">\n",
       "    <rect x=\"384\" y=\"0\" width=\"1345\" height=\"1280\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip3700)\" d=\"\n",
       "M137.065 1108.82 L1872.76 1108.82 L1872.76 121.675 L137.065 121.675  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip3702\">\n",
       "    <rect x=\"137\" y=\"121\" width=\"1737\" height=\"988\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip3702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  368.741,1108.82 368.741,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  733.847,1108.82 733.847,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1098.95,1108.82 1098.95,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1464.06,1108.82 1464.06,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1829.16,1108.82 1829.16,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  137.065,986.016 1872.76,986.016 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  137.065,837.791 1872.76,837.791 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  137.065,689.566 1872.76,689.566 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  137.065,541.341 1872.76,541.341 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  137.065,393.116 1872.76,393.116 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  137.065,244.891 1872.76,244.891 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  137.065,1108.82 1872.76,1108.82 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  137.065,1108.82 137.065,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  368.741,1108.82 368.741,1094.01 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  733.847,1108.82 733.847,1094.01 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1098.95,1108.82 1098.95,1094.01 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1464.06,1108.82 1464.06,1094.01 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1829.16,1108.82 1829.16,1094.01 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  137.065,986.016 163.1,986.016 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  137.065,837.791 163.1,837.791 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  137.065,689.566 163.1,689.566 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  137.065,541.341 163.1,541.341 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  137.065,393.116 163.1,393.116 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  137.065,244.891 163.1,244.891 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 368.741, 1160.42)\" x=\"368.741\" y=\"1160.42\">-2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 733.847, 1160.42)\" x=\"733.847\" y=\"1160.42\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1098.95, 1160.42)\" x=\"1098.95\" y=\"1160.42\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1464.06, 1160.42)\" x=\"1464.06\" y=\"1160.42\">4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1829.16, 1160.42)\" x=\"1829.16\" y=\"1160.42\">6</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 1003.52)\" x=\"117.865\" y=\"1003.52\">-5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 855.291)\" x=\"117.865\" y=\"855.291\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 707.066)\" x=\"117.865\" y=\"707.066\">5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 558.841)\" x=\"117.865\" y=\"558.841\">10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 410.616)\" x=\"117.865\" y=\"410.616\">15</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 262.391)\" x=\"117.865\" y=\"262.391\">20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:84px; text-anchor:middle;\" transform=\"rotate(0, 1004.91, 73.2)\" x=\"1004.91\" y=\"73.2\">fn1() and its derivative</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1004.91, 1237.88)\" x=\"1004.91\" y=\"1237.88\">x</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip3702)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  186.188,354.281 197.252,368.905 208.316,383.311 219.379,397.499 230.443,411.469 241.507,425.222 252.571,438.757 263.635,452.074 274.698,465.173 285.762,478.055 \n",
       "  296.826,490.719 307.89,503.165 318.954,515.393 330.017,527.404 341.081,539.196 352.145,550.771 363.209,562.128 374.273,573.268 385.337,584.189 396.4,594.893 \n",
       "  407.464,605.379 418.528,615.647 429.592,625.698 440.656,635.53 451.719,645.145 462.783,654.542 473.847,663.722 484.911,672.683 495.975,681.427 507.039,689.953 \n",
       "  518.102,698.261 529.166,706.352 540.23,714.224 551.294,721.879 562.358,729.316 573.421,736.536 584.485,743.537 595.549,750.321 606.613,756.887 617.677,763.235 \n",
       "  628.741,769.366 639.804,775.278 650.868,780.973 661.932,786.45 672.996,791.71 684.06,796.751 695.123,801.575 706.187,806.181 717.251,810.569 728.315,814.74 \n",
       "  739.379,818.692 750.442,822.427 761.506,825.945 772.57,829.244 783.634,832.325 794.698,835.189 805.762,837.835 816.825,840.263 827.889,842.474 838.953,844.467 \n",
       "  850.017,846.241 861.081,847.799 872.144,849.138 883.208,850.259 894.272,851.163 905.336,851.849 916.4,852.317 927.464,852.568 938.527,852.601 949.591,852.415 \n",
       "  960.655,852.013 971.719,851.392 982.783,850.553 993.846,849.497 1004.91,848.223 1015.97,846.731 1027.04,845.022 1038.1,843.095 1049.17,840.949 1060.23,838.587 \n",
       "  1071.29,836.006 1082.36,833.207 1093.42,830.191 1104.48,826.957 1115.55,823.505 1126.61,819.836 1137.68,815.948 1148.74,811.843 1159.8,807.52 1170.87,802.98 \n",
       "  1181.93,798.221 1193,793.245 1204.06,788.051 1215.12,782.639 1226.19,777.01 1237.25,771.162 1248.31,765.097 1259.38,758.814 1270.44,752.314 1281.51,745.595 \n",
       "  1292.57,738.659 1303.63,731.505 1314.7,724.133 1325.76,716.544 1336.82,708.736 1347.89,700.711 1358.95,692.468 1370.02,684.008 1381.08,675.329 1392.14,666.433 \n",
       "  1403.21,657.319 1414.27,647.987 1425.34,638.438 1436.4,628.67 1447.46,618.685 1458.53,608.482 1469.59,598.062 1480.65,587.423 1491.72,576.567 1502.78,565.493 \n",
       "  1513.85,554.201 1524.91,542.692 1535.97,530.964 1547.04,519.019 1558.1,506.856 1569.16,494.476 1580.23,481.877 1591.29,469.061 1602.36,456.027 1613.42,442.775 \n",
       "  1624.48,429.305 1635.55,415.618 1646.61,401.713 1657.68,387.59 1668.74,373.249 1679.8,358.691 1690.87,343.915 1701.93,328.921 1712.99,313.709 1724.06,298.279 \n",
       "  1735.12,282.632 1746.19,266.767 1757.25,250.684 1768.31,234.383 1779.38,217.865 1790.44,201.129 1801.5,184.175 1812.57,167.003 1823.63,149.613 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3702)\" style=\"stroke:#ff0000; stroke-width:4; stroke-opacity:1; fill:none\" stroke-dasharray=\"2, 4\" points=\"\n",
       "  186.188,1080.88 197.252,1077.29 208.316,1073.69 219.379,1070.1 230.443,1066.51 241.507,1062.91 252.571,1059.32 263.635,1055.73 274.698,1052.13 285.762,1048.54 \n",
       "  296.826,1044.95 307.89,1041.35 318.954,1037.76 330.017,1034.17 341.081,1030.57 352.145,1026.98 363.209,1023.39 374.273,1019.79 385.337,1016.2 396.4,1012.61 \n",
       "  407.464,1009.01 418.528,1005.42 429.592,1001.83 440.656,998.234 451.719,994.641 462.783,991.047 473.847,987.454 484.911,983.86 495.975,980.267 507.039,976.674 \n",
       "  518.102,973.08 529.166,969.487 540.23,965.894 551.294,962.3 562.358,958.707 573.421,955.114 584.485,951.52 595.549,947.927 606.613,944.334 617.677,940.74 \n",
       "  628.741,937.147 639.804,933.554 650.868,929.96 661.932,926.367 672.996,922.774 684.06,919.18 695.123,915.587 706.187,911.994 717.251,908.4 728.315,904.807 \n",
       "  739.379,901.214 750.442,897.62 761.506,894.027 772.57,890.434 783.634,886.84 794.698,883.247 805.762,879.654 816.825,876.06 827.889,872.467 838.953,868.874 \n",
       "  850.017,865.28 861.081,861.687 872.144,858.094 883.208,854.5 894.272,850.907 905.336,847.314 916.4,843.72 927.464,840.127 938.527,836.534 949.591,832.94 \n",
       "  960.655,829.347 971.719,825.754 982.783,822.16 993.846,818.567 1004.91,814.974 1015.97,811.38 1027.04,807.787 1038.1,804.194 1049.17,800.6 1060.23,797.007 \n",
       "  1071.29,793.414 1082.36,789.82 1093.42,786.227 1104.48,782.634 1115.55,779.04 1126.61,775.447 1137.68,771.854 1148.74,768.26 1159.8,764.667 1170.87,761.074 \n",
       "  1181.93,757.48 1193,753.887 1204.06,750.294 1215.12,746.7 1226.19,743.107 1237.25,739.514 1248.31,735.92 1259.38,732.327 1270.44,728.734 1281.51,725.14 \n",
       "  1292.57,721.547 1303.63,717.954 1314.7,714.36 1325.76,710.767 1336.82,707.174 1347.89,703.58 1358.95,699.987 1370.02,696.394 1381.08,692.8 1392.14,689.207 \n",
       "  1403.21,685.614 1414.27,682.02 1425.34,678.427 1436.4,674.834 1447.46,671.24 1458.53,667.647 1469.59,664.054 1480.65,660.46 1491.72,656.867 1502.78,653.273 \n",
       "  1513.85,649.68 1524.91,646.087 1535.97,642.493 1547.04,638.9 1558.1,635.307 1569.16,631.713 1580.23,628.12 1591.29,624.527 1602.36,620.933 1613.42,617.34 \n",
       "  1624.48,613.747 1635.55,610.153 1646.61,606.56 1657.68,602.967 1668.74,599.373 1679.8,595.78 1690.87,592.187 1701.93,588.593 1712.99,585 1724.06,581.407 \n",
       "  1735.12,577.813 1746.19,574.22 1757.25,570.627 1768.31,567.033 1779.38,563.44 1790.44,559.847 1801.5,556.253 1812.57,552.66 1823.63,549.067 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3702)\" style=\"stroke:#0000ff; stroke-width:4; stroke-opacity:1; fill:none\" stroke-dasharray=\"16, 10\" points=\"\n",
       "  186.188,1080.88 197.252,1077.29 208.316,1073.69 219.379,1070.1 230.443,1066.51 241.507,1062.91 252.571,1059.32 263.635,1055.73 274.698,1052.13 285.762,1048.54 \n",
       "  296.826,1044.95 307.89,1041.35 318.954,1037.76 330.017,1034.17 341.081,1030.57 352.145,1026.98 363.209,1023.39 374.273,1019.79 385.337,1016.2 396.4,1012.61 \n",
       "  407.464,1009.01 418.528,1005.42 429.592,1001.83 440.656,998.234 451.719,994.641 462.783,991.047 473.847,987.454 484.911,983.86 495.975,980.267 507.039,976.674 \n",
       "  518.102,973.08 529.166,969.487 540.23,965.894 551.294,962.3 562.358,958.707 573.421,955.114 584.485,951.52 595.549,947.927 606.613,944.334 617.677,940.74 \n",
       "  628.741,937.147 639.804,933.554 650.868,929.96 661.932,926.367 672.996,922.774 684.06,919.18 695.123,915.587 706.187,911.994 717.251,908.4 728.315,904.807 \n",
       "  739.379,901.214 750.442,897.62 761.506,894.027 772.57,890.434 783.634,886.84 794.698,883.247 805.762,879.654 816.825,876.06 827.889,872.467 838.953,868.874 \n",
       "  850.017,865.28 861.081,861.687 872.144,858.094 883.208,854.5 894.272,850.907 905.336,847.314 916.4,843.72 927.464,840.127 938.527,836.534 949.591,832.94 \n",
       "  960.655,829.347 971.719,825.754 982.783,822.16 993.846,818.567 1004.91,814.974 1015.97,811.38 1027.04,807.787 1038.1,804.194 1049.17,800.6 1060.23,797.007 \n",
       "  1071.29,793.414 1082.36,789.82 1093.42,786.227 1104.48,782.634 1115.55,779.04 1126.61,775.447 1137.68,771.854 1148.74,768.26 1159.8,764.667 1170.87,761.074 \n",
       "  1181.93,757.48 1193,753.887 1204.06,750.294 1215.12,746.7 1226.19,743.107 1237.25,739.514 1248.31,735.92 1259.38,732.327 1270.44,728.734 1281.51,725.14 \n",
       "  1292.57,721.547 1303.63,717.954 1314.7,714.36 1325.76,710.767 1336.82,707.174 1347.89,703.58 1358.95,699.987 1370.02,696.394 1381.08,692.8 1392.14,689.207 \n",
       "  1403.21,685.614 1414.27,682.02 1425.34,678.427 1436.4,674.834 1447.46,671.24 1458.53,667.647 1469.59,664.054 1480.65,660.46 1491.72,656.867 1502.78,653.273 \n",
       "  1513.85,649.68 1524.91,646.087 1535.97,642.493 1547.04,638.9 1558.1,635.307 1569.16,631.713 1580.23,628.12 1591.29,624.527 1602.36,620.933 1613.42,617.34 \n",
       "  1624.48,613.747 1635.55,610.153 1646.61,606.56 1657.68,602.967 1668.74,599.373 1679.8,595.78 1690.87,592.187 1701.93,588.593 1712.99,585 1724.06,581.407 \n",
       "  1735.12,577.813 1746.19,574.22 1757.25,570.627 1768.31,567.033 1779.38,563.44 1790.44,559.847 1801.5,556.253 1812.57,552.66 1823.63,549.067 \n",
       "  \"/>\n",
       "<path clip-path=\"url(#clip3700)\" d=\"\n",
       "M614.515 418.315 L1159.58 418.315 L1159.58 176.395 L614.515 176.395  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  614.515,418.315 1159.58,418.315 1159.58,176.395 614.515,176.395 614.515,418.315 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  633.715,236.875 748.915,236.875 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 768.115, 254.375)\" x=\"768.115\" y=\"254.375\">fn1()</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#ff0000; stroke-width:4; stroke-opacity:1; fill:none\" stroke-dasharray=\"2, 4\" points=\"\n",
       "  633.715,297.355 748.915,297.355 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 768.115, 314.855)\" x=\"768.115\" y=\"314.855\">crude derivative</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip3700)\" style=\"stroke:#0000ff; stroke-width:4; stroke-opacity:1; fill:none\" stroke-dasharray=\"16, 10\" points=\"\n",
       "  633.715,357.835 748.915,357.835 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip3700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 768.115, 375.335)\" x=\"768.115\" y=\"375.335\">ForwardDiff pkg </text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot( [x x x],[fn1.(x,0.5) dydx_A dydx_B],\n",
    "      line = [:solid :dot :dash],\n",
    "      linecolor = [:black :red :blue],\n",
    "      label = [\"fn1()\" \"crude derivative\" \"ForwardDiff pkg \"],\n",
    "      legend = :top,\n",
    "      title = \"fn1() and its derivative\",\n",
    "      xlabel = \"x\",\n",
    "      ylabel = \"\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Comments on Using the ForwardDiff Package (extra)\n",
    "\n",
    "The `ForwardDiff` package applies an interesting approach to calculate derivatives, using a special number type (\"dual numbers\"). This means that your code must be able to handle such numbers. In most cases, that is not a problem, but you may have to watch out if you create arrays to store (intermediate?) results inside the function. See the examples below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "fnDoesWork (generic function with 1 method)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function fnDoesNotWork(b,a)\n",
    "    z = zeros(length(b))              #will not work with ForwardDiff, since\n",
    "    for i = 1:length(z)               #z cannot store dual numbers\n",
    "        z[i] = b[i]*i\n",
    "    end\n",
    "    return sum(z) + a\n",
    "end\n",
    "\n",
    "function fnDoesWork(b,a)\n",
    "    z = zeros(eltype(b),length(b))   #will work with ForwardDiff, since\n",
    "    for i = 1:length(z)              #when b is a dual number, so is z \n",
    "        z[i] = b[i]*i                #could also start with z = similar(b)\n",
    "    end\n",
    "    return sum(z) + a\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 1.0\n",
       " 2.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b0 = [1.5,2]\n",
    "\n",
    "#ForwardDiff.gradient(b->fnDoesNotWork(b,1),b0)      #uncomment to get an error\n",
    "ForwardDiff.gradient(b->fnDoesWork(b,1),b0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 1.3.0",
   "language": "julia",
   "name": "julia-1.3"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.3.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
